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Abstract 

We demonstrate a possibility to stabilize three-dimensional spatiotemporal solitons 
("light bullets") in self-focusing Kerr media by means of a combination of disper- 
sion management in the longitudinal direction (with the group-velocity dispersion 
alternating between positive and negative values) and periodic modulation of the 
refractive index in one transverse direction, out of the two. The analysis is based 
on the variational approximation (results of direct three-dimensional simulations 
will be reported in a follow-up work). A predicted stability area is identified in the 
model's parameter space. It features a minimum of the necessary strength of the 
transverse modulation of the refractive index, and finite minimum and maximum 
values of the soliton's energy. The former feature is also explained analytically. 



1 Introduction 



Search for spatiotemporal solitons in diverse optical media, alias "light bullets" 
(LBs) [1], is a challenge to fundamental and applied research in nonlinear 
optics, see original works [2,3,4,5,6,7,8,9,10] and a very recent review [11]. 
Stationary solutions for LBs can easily be found in the cubic (x ) multi- 
dimensional nonlinear Schrodinger (NLS) equation [1], but their stability is 
a problem, as they are unstable against spatiotemporal collapse [12]. The 
problem may be avoided by resorting to milder nonlinearities, such as saturable 
[9], cubic-quintic [13], or quadratic (x^) [2,3,4,5,6,7,8]. 
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Despite considerable progress in theoretical studies, three-dimensional (3D) 
LBs in a bulk medium have not yet been observed in an experiment. The 
only successful experimental finding reported thus far was a stable quasi-2D 
spatiotemporal soliton in x^ crystals [7] (the tilted- wavefront technique [14], 
used in that work, precluded achieving self- confinement in one transverse di- 
rection). On the other hand, it was predicted [8] that a spatial cylindrical 
soliton may be stabilized in a bulk medium composed of layers with alternat- 
ing signs of the Kerr coefficient. Similar stabilization was then predicted for 
what may be regarded as 2D solitons in Bose-Einstein condensates (BECs), 
with the coefficient in front of the cubic nonlinear term subjected to periodic 
modulation in time via the Feshbach resonance in external ac magnetic field 
[15,16]. However, no stable 3D soliton could be found in either realization 
(optical or BEC) of this setting. 

Serious difficulties encountered in the experimental search for LBs in 3D media 
is an incentive to look for alternative settings admitting stable 3D optical soli- 
tons. With the Kerr nonlinearity, a possibility is to use a layered structure that 
periodically reverses the sign of the local group- velocity dispersion (GVD), 
without affecting the x^ coefficient. This resembles a well-known scheme in 
fiber optics, known as dispersion management (DM), see, e.g., Refs. [17] and 
review [18]. A 2D generalization of the DM scheme was recently proposed, 
assuming a layered planar waveguide of this type, uniform in the transverse 
direction [19,20]. As a result, large stability regions for the 2D spatiotempo- 
ral solitons were identified, including double-peaked breathers; however, a 3D 
version of the same model could not give rise to any stable soliton [19]. It was 
also shown in Ref. [19] that no stable 3D soliton could be found in a more 
sophisticated model, which combines the DM and periodic modulation of the 
Kerr coefficient in the longitudinal direction. 

Another approach to the stabilization of multidimensional solitons was devel- 
oped in the context of the self- attracting BEC. It is based on the correspond- 
ing Gross-Pitaevskii equation which includes a periodic potential created as 
an optical lattice (OL, i.e., an interference pattern produced by illuminating 
the condensate by counter-propagating coherent laser beams). It has been 
demonstrated that 2D [21,22,23] and 3D [21] solitons can be easily stabilized 
by the OL of the same dimension. Moreover, stable solitons can also be readily 
supported by low- dimensional OLs, i.e., ID and 2D ones in the 2D [23,24] and 
3D [23,24,25] cases, respectively; additionally, a 3D soliton can be stabilized 
by a cylindrical (Bessel) lattice [26], similar to one introduced, in the context 
of 2D models, in Ref. [27] . On the other hand, 3D solitons cannot be stabilized 
by a ID periodic potential [24]; however, the ID lattice potential in combina- 
tion with the above-mentioned time-periodic modulation of the nonlinearity, 
provided by the Feshbach resonance in the ac magnetic field, supports single- 
and multi-peaked stable 3D solitons in vast areas of the respective parameter 
space [28]. 
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The above results suggest a possibility of existence of stable 3D "bullets" in 
a medium with the DM in the longitudinal direction (z), additionally 
equipped with an effective lattice potential (i.e., periodic modulation of the 
refractive index) in one transverse direction (y), while in the remaining trans- 
verse direction (x) the medium remains uniform. If this is possible, stable LBs 
will be definitely possible too in a medium with the periodic modulation of 
the refractive index in both transverse directions; however, the setting with 
one uniform direction is more interesting in terms of steering solitons and 
studying collisions between them [23,24]. The objective of the present work is 
to predict such 3D spatiotemporal solitons and investigate their stability. Our 
first consideration of this possibility is based on the variational approximation 
(VA); systematic simulations of the 3D model are quite complicated, and will 
be presented in a follow-up work. It is relevant to mention that the existence 
and stability of 3D solitons in the Gross-Pitaevskii equation with the quasi- 
2D periodic potential, which were originally predicted by the VA[23,24], was 
definitely confirmed by direct simulations [23,24,25], which suggests that in 
the present model the 3D solitons may easily be stable too. 

The model is based on the normalized NLS equation describing the evolution of 
the local amplitude u of the electromagnetic wave, which is a straightforward 
extension of the 2D model put forward in Ref. [19]: 



.du 

dz 



' d 2 d 2 d 2 \ 

w + W + D{z) ^) +£CO<2y) + H 



u = 0. (1) 



Here, e is the strength of the transverse modulation (the modulation period is 
normalized to be 7r), while r and D(z) are the same reduced temporal variable 
and local GVD coefficient as in the fiber-optic DM models [17,18]. 

Equation (1) implies the propagation of a linearly polarized wave, with the 
single component u; a more general situation will be described by a two- 
component (vectorial) version of Eq. (1), with the two polarization coupled, 
as usual, by the cubic cross-phase-modulation terms. We do not expect that 
the vectorial model will produce results qualitatively different form those pre- 
sented below. As usual, the NLS equation assumes the applicability of the 
paraxial approximation, i.e., the spatial size of solitons (see below) must be 
much larger than the underlying wavelength of light, which is definitely a 
physically relevant assumption [11], and the temporal part of the equation 
implies that the higher-order GVD is negligible (previous considerations have 
demonstrated that the higher-order dispersion does not drastically alter DM 
solitons [29]). 

As is commonly adopted, we assume a symmetric DM map, with equal lengths 
L of the normal- and anomalous-GVD segments (usually, the results are not 
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sensitive to the map's asymmetry), 

D + D m > 0, < z < L, 



D(z) = 



(2) 

D - An < 0, L < z < 2L, 



;he average dispersion being much smaller than the modulation amplitude, 
D <C An- Using the scaling invariances of Eq. (1), we fix L = 1 and D m = I. 

Recently, a somewhat similar 2D model was introduced in Ref. [30]. The most 
important difference is that it has the variable coefficient D(z) multiplying 
both the GVD and diffraction terms, u TT and u xx . Actually, that model was 
motivated by a continuum limit of some discrete systems; in the present con- 
text, it would be quite difficult to implement the periodic reversal of the sign 
of the transverse diffraction. 



2 The variational approximation 



Aiming to apply the VA for the search of LB solutions (a review of the varia- 
tional method can be found in Ref. [18]), we adopt the Gaussian ansatz, 



u = A(z) exp I i(f)(z) - - 



x 



+ 



+ 



W 2 (z) V 2 (z) T 2 (z) 



+ 



+ 2 



b(z) x 2 + c(z)y 2 + P(z) r 2 



(3) 



where A and are the amplitude and phase of the soliton, T and W, V are its 
temporal and two transverse spatial widths, and (3 and b, c are the temporal 
and two spatial chirps. The Lagrangian from which Eq. (1) can be derived is 
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oo +oo +oo 

dx J dy J dr 

— oo — oo — oo 

+2ecos(2?/)H 2 + H 4 



*/ * *\ I 1 2 i 1 2 T-^ii 1 2 

i [u z u — u z u) — \u x \ — \Uy\ — D \u T \ 



(4) 
(5) 



The substitution of the ansatz (3) in this expression and integrations lead to 
an effective Lagrangian, with the prime standing for d/dz: 



(4/tt 3/2 ) L cS = A 2 WVT [40' - b'W 2 - c'V 2 - f3'T 2 - W~ 2 - V~ 2 - DT~ 2 
-b 2 W 2 - c 2 V 2 + e exp (-V 2 ) - D(z)p 2 T 2 + A 2 /V2] , (6) 
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The first variational equation, 5L c s/5(p = 0, applied to Eq. (6) yields the 
energy conservation, dE/dz = 0, with 

E = A 2 WVT. (7) 



The conservation of E is used to eliminate A 2 from the set of subsequent 
equations, SL eS /S (W, V, T, b, c, (3) = 0. They can be arranged so as, first, to 
eliminate the chirps, 

b = W'/W, c = V'/V, (3 = D^T'/T. (8) 



the remaining equations for the spatial and temporal widths being 



1 E 

W" = — - (9) 

V"=-^-4eVexp(-V 2 ) =^ , (10) 

V 3 P V y 2y/2WV 2 T 

'T'Y _ D E 
D ~T3 " 2 ^2WVT 2 ' 
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The Hamiltonian of these equations, which is a dynamical invariant in the 
case of constant D, is 

^ , ^ , (Tf , 1 , 1 , D A _, _ E 



n = (wy + tyy + ^-f - + — + _ + — -45 ex P (-\/ 2 ) - 



D W 2 V 2 T 2 > ^WVT 



In the case of the piece-wise constant modulation, such as in (2), the variables 
W, W, V, V, T and (3 must be continuous at junctions between the segments 
with D± = D ±D m . As it follows from Eq. (8), the continuity of the temporal 
chirp /3(z) implies a jump of the derivative T' when passing from D_ to D + , 
or vice versa: 

(T') D=D+ = (D + /D_) (T') D=D _ . (12) 



In the case of a continuous DM map, rather than the one (2), Eq. (11) has 
a formal singularity at the points where D(z) vanishes, changing its sign. 
However, it is known that there is no real singularity in this T' vanishes 

at the same points, which cancels the singularity out [18]. 

In the absence of the DM and transverse modulation, i.e., D = +1 and e — 0, 
three equations (9) - (11) reduce to one, which is tantamount to the variational 
equation derived in Ref. [31] from the spatiotemporally isotropic ansatz [cf. 
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Eq. (3)], u = A exp [i0 - (1/2) (W~ 2 + ife) (x 2 + y 2 + r 2 )]. In particular, this 
single equation correctly predicts the asymptotic law of the strong collapse in 
the 3D case, which is stable against anisotropic perturbations [32], V = W = 

T « (5E/3v^) (zo — -z) 2 ^ 5 , z = zq being the collapse point. The location of 
this point is determined by initial conditions, but, in any case, it belongs to 
an interval D > 0, where the GVD is anomalous. 

Another possible collapse scenario is an effectively two-dimensional (weak) 
one, with two widths shrinking to zero as z — z — > 0, while the third one 
remains finite. For instance, the corresponding asymptotic law may be 



V ; 4 + A 4 4V2EA 2 



[z - z) In (zq - z) , (13) 



where A is a positive constant or else V = 2/T ~ (z — z) 1 ^ 2 and W — > W 
(in this case too, the collapse point Zq must belong to a segment with D > 0). 
In direct simulations of Eqs. (9) - (11), we actually observed only the latter 
scenario. However, we did not specially try to find initial conditions that could 
initiate a solution corresponding to the strong 3D collapse, as our objective 
is not the study of the collapse, but rather search for solitons stable against 
collapse. In fact, known results for the solitons in the 3D Gross-Pitaevskii 
equation with the OL potential suggest that, while the VA may be incorrect 
in the description of the collapse, as a singular solution, it provides for quite 
accurate predictions for the stability of solitons as regular solutions [23,24]. 

If the DM is absent, and the constant GVD is normal, i.e., D = — 1, only the 
2D collapse in the transverse plane would be possible, so that (cf. Eq. (13)) 

V = W = A (z - z) l/2 , T = -^L + ( A + Ai ) , _ n ln ( _ ) _ 
v ; 4 + A 4 4V2EA 2 V ; V ; 

However, we did not observed this collapse scenario in our simulations. The 
same comment as one given above pertains to this case as well. 

A possibility of the stabilization of the 3D soliton by a sufficiently strong 
lattice can be understood noticing that, for large e, one may keep only the 
first two terms on the right-hand side of Eq. (10). This approximation yields 
a nearly constant value Vq of V, which is a smaller root of the corresponding 
equation, 

4eV 4 exp (-V 2 ) = 1 (14) 
(a larger root corresponds to an unstable equilibrium). The two roots exist 
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provided that 



e > e- 



mm 



e 2 /16 w 0.46, 



(15) 



the relevant one being limited by V < 2. Then, the substitution of V — V 
in the remaining equations (9) and (11) leads to essentially the same VA- 
generated dynamical system as derived for the 2D DM model in Ref. [19], 
which was shown to give rise to stable spatiotemporal solitons. On the other 
hand, it was demonstrated in Ref. [19] too that, in the case of e = 0, the 3D 
VA equations, as well as the full underlying 3D model, have no stable soliton 
solutions. 

The stabilization of the LB in the present model for large e can also be under- 
stood in a different way, without resorting to VA: in a very strong lattice, the 
soliton is trapped entirely in a single "valley" of the periodic potential, and 
the problem thus reduces to a nearly 2D one, where spatiotemporal solitons 
may be stable, cf. a similar stabilization mechanism for the solitons in the 
Gross-Pitaevskii equations developed in [16]. From this point of view, a really 
interesting issue is to find an actual minimum e min of the lattice's strength 
which is necessary for the stabilization of the 3D solitons, as at e close enough 
to £ min the stabilized solitons are truly 3D objects, rather than their nearly-2D 
counterparts. 



3 Results 



We explored the parameter space of the variational system (9) - (11), [E,e, Dj, 

by means of direct simulations of the equations (with regard to the jump con- 
dition (12)). As a result, it was possible to identify regions where the model 
admits stable solitons featuring regular oscillations in z with the DM-map 
period. An example of such a regime is shown in Fig. 1 (oscillations in the 
evolution of W are not visible in the figure because, as an estimate demon- 
strates, their amplitude is ~ 0.001). 

Systematic results obtained from the simulations are summarized in stability 
diagrams displayed in Figs. 2 and 3. A remarkable fact, apparent in Fig. 2, 
is that the minimum value of the lattice's strength, e min = 0.46, at which the 
solitons may be stable, coincides with the analytical prediction (15), up to the 
available numerical accuracy. 

The existence of a maximum value E max of the energy admitting the stable LBs 
is, essentially, a quasi-2D feature, which can be understood assuming that the 
potential lattice is strong. Indeed, as explained above, in such a case the value 
of V is approximately fixed as the smaller root of Eq. (14). Within a segment 
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Fig. 1. An example of the stable evolution of solutions to the variational equations 
(9) - (11). The soliton's widths in the direction x,y and r, i.e., W, V and T, are 
shown as functions of z, for E = 0.5, e = 1, and D = 0. 



where the GVD coefficient keeps the constant value, D + = D + D m > 0, which 
corresponds to anomalous dispersion (see Eq. (2), the remaining equations (9) 
and (11) are tantamount to those for a uniform 2D Kerr-self-focusing medium, 
hence the energy is limited by the value -E max corresponding to the Townes 
soliton; the soliton will collapse if E > -E max [12]. 



The fact that the region of stable solitons is also limited by a minimum energy, 
E min , except for the case of D = 0, when E min = (see Fig. 3)), is actually a 
quasi- ID feature, which is characteristic to the DM solitons in optical fibers. 
In that case, the term ~ E in the evolution equation for T(z), cf. Eq. (11), 
is necessary to balance the average GVD coefficient D, so that E min and D 
vanish simultaneously [17]. It is noteworthy too that, as well as in the case 
of the ID DM solitons in fibers, the stability area in Fig. 3 includes a part 
with normal average GVD, D < 0, which seems counterintuitive, but can be 
explained [17]. This part extends in Fig. 3 up to \ —D) m 0.005. 

V / max 
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1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

0.5 1 1.5 

E 

Fig. 2. The stability area for the 3D spatiotemporal solitons in the (E, e) plane, 
with D = 0, is shown by light-gray shading. In gray and dark-gray regions, the 
3D soliton is predicted, respectively, to spread out and collapse. The vertical line 
corresponds to the analytically predicted threshold (15). 

4 Conclusions 



In this work, we have proposed a possibility to stabilize spatiotemporal solitons 
( "light bullets" ) in three-dimensional self-focusing Kerr media by means of the 
dispersion management (DM), which means that the local group- velocity dis- 
persion coefficient alternates between positive and negative values along the 
propagation direction, z. Recently, it was shown that the DM alone can stabi- 
lize solitons in 2D (planar) waveguides, but in the bulk (3D) DM medium the 
"bullets" are unstable. In this work, we have demonstrated that the complete 
stabilization can be provided if the longitudinal DM is combined with peri- 
odic modulation of the refractive index in one transverse direction (y), out of 
the two. The analysis was based on the variational approximation (systematic 
results of direct simulations will be reported in a follow-up paper). A stability 
area for the light bullets was identified in the model's parameter space. Its 
salient features are a necessary minimum strength of the transverse modu- 
lation of the refractive index, and minimum and maximum values E min ^ max 
of the soliton's energy. The former feature can be accurately predicted (see 
Eq. (15)) in an analytical form from the evolution equation for the width of 
the soliton in the y-direction. The existence of E min , which vanishes when we 
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Fig. 3. The stability area in the (E,D) pi ane, with e = 1, is shown by light-gray 
shading. In the gray region, the 3D soliton is predicted to spread out. 

assume zero average dispersion, can be explained in the same way as for the 
temporal solitons in DM optical fibers. Also, similar to the case of DM solitons 
in fibers, we find that the stability area extends to a region of normal average 
dispersion [17]. On the other hand, the existence of E max can be understood 
similarly to as it was recently done in the 2D counterpart of the present model 
(the strong transverse lattice can squeeze the system to a nearly 2D shape). 

The results presented in this work suggest a new approach to the challenging 
problem of the creation of 3D spatiotemporal optical solitons. The model 
also opens a way to address advanced issues, such as collisions between the 
LBs, and the existence and stability of solitons with different symmetries (for 
instance, LBs which are odd in the longitudinal and/or transverse directions). 
These issues will be considered elsewhere. 
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